COMMENT ⊗   VALID 00003 PAGES 
RECORD PAGE   DESCRIPTION
 00001 00001
 00002 00002	ENTRY DECOM,SOLVE,QSOLVE,INVRT,DETERM
 00006 00003	SIMPLE INTERNAL PROCEDURE SOLVE(REAL ARRAY LUA,B,XINTEGER ARRAY PS)
 00010 ENDMK
⊗;
ENTRY DECOM,SOLVE,QSOLVE,INVRT,DETERM;
BEGIN "MATRIX" COMMENT  This is the Forsythe-Moler matrix package (op.cit. sec 16).
	IMPROVE (see pp60-61) can be added later if more accuracy is desired;
DEFINE ⊃="COMMENT",STX="1 STEP 1 UNTIL",CRLF="'15&'12",
	ERROR(X)="BEGIN OUTSTR(CRLF&""X"");RETURN END";

SIMPLE PROCEDURE SINGULAR(INTEGER WHY); ⊃ Prints error messages;
	BEGIN IF WHY=0 THEN OUTSTR("MATRIX WITH ZERO ROW IN DECOM")
		ELSE IF WHY=1 THEN OUTSTR("SINGULAR MATRIX IN DECOM.
			SOLVE WILL DIVIDE BY ZERO");
	END "SINGULAR";

INTERNAL PROCEDURE DECOM(REAL ARRAY A,LUA;INTEGER ARRAY PS);
		⊃ This is F-M DECOMPOSE. It performs Gauss-Jordan
		reduction with scaling and interchanging of rows
		yielding an intermediate matrix LUA. LUA is the 
		matrix whose lower triangle L with a main diagonal
		of 1's, and whose upper triangle together with the
		main diagonal satisfy LU=PA. L contains the multipliers
		M[i,j]=A[i,j]/A[j,j]. PS is the permuation table
		from which a permutation matrix can be constructed
		such that LU=PA;
	BEGIN INTEGER I,J,K,PIVIND,N; REAL NRMROW,PIVOT,SIZE,BIGGST,MULT,TEM;
	N←ARRINFO(A,2); IF ARRINFO(A,4)≠N THEN ERROR(DECOM: [A] NOT SQUARE);
	IF ¬(ARRINFO(LUA,2)=ARRINFO(LUA,4)=N) THEN ERROR(DECOM: LUA WRONG SIZE);
	IF ARRINFO(PS,2)<N THEN ERROR(DECOM: PS TOO SMALL);
	BEGIN "COMPUT" REAL ARRAY SCALES[1:N]; LABEL ENDKLOOP;
⊃  Initialize LUA,PS,SCALES; ARRTRAN(LUA,A);
	FOR I←STX N DO BEGIN PS[I]←I;NRMROW←0;
		FOR J←STX N DO IF NRMROW<(TEM←ABS(LUA[I,J])) THEN NRMROW←TEM;
		IF NRMROW≠0 THEN SCALES[I]←1/NRMROW
		ELSE BEGIN SCALES[I]←0;SINGULAR(0) END;END;
⊃ Guassian elimination with partial pivoting;
	FOR K←STX N-1 DO BEGIN BIGGST←0;
		FOR I←K STEP 1 UNTIL N DO BEGIN
			SIZE←ABS(LUA[PS[I],K])*SCALES[PS[I]];
			IF BIGGST<SIZE THEN BEGIN BIGGST←SIZE;PIVIND←I END;END;
		IF BIGGST=0 THEN BEGIN SINGULAR(1);GO ENDKLOOP END;
		IF PIVIND≠K THEN PS[K]↔PS[PIVIND];
		PIVOT←LUA[PS[K],K];
		FOR I←K+1 STEP 1 UNTIL N DO BEGIN
			MULT←LUA[PS[I],K]←LUA[PS[I],K]/PIVOT;
			IF MULT≠0 THEN FOR J←K+1 STEP 1 UNTIL N DO
				LUA[PS[I],J]←LUA[PS[I],J]-MULT*LUA[PS[K],J];END;
	ENDKLOOP: END;
	IF LUA[PS[N],N]=0 THEN SINGULAR(1);
	END "COMPUT"; END "DECOM";
SIMPLE INTERNAL PROCEDURE SOLVE(REAL ARRAY LUA,B,X;INTEGER ARRAY PS);
		⊃ Solves AX=B using LUA and PS generated by DECOM;
	BEGIN INTEGER I,J,N; REAL DOT; N←ARRINFO(LUA,2);
	IF ARRINFO(LUA,4)≠N THEN ERROR(SOLVE: [LUA] NOT SQUARE);
	IF ARRINFO(B,2)<N THEN ERROR(SOLVE: B TOO SMALL);
	IF ARRINFO(X,2)<N THEN ERROR(SOLVE: X TOO SMALL);
	FOR I←STX N DO BEGIN DOT←0;
		FOR J←STX I-1 DO DOT←DOT+LUA[PS[I],J]*X[J];
		X[I]←B[PS[I]]-DOT END;
	FOR I←N STEP -1 UNTIL 1 DO BEGIN DOT←0;
		FOR J←I+1 STEP 1 UNTIL N DO DOT←DOT+LUA[PS[I],J]*X[J];
		X[I]←(X[I]-DOT)/LUA[PS[I],I] END;
	END "SOLVE";

INTERNAL PROCEDURE QSOLVE(REAL ARRAY A,X,B);
		⊃ Quick solution to AX=B. Does not save LUA and PS.
		Calls DECOM and SOLVE;
	BEGIN INTEGER N; N←ARRINFO(A,2);
	IF ARRINFO(A,4)≠N THEN ERROR(QSOLVE: [A] NOT SQUARE);
	IF ARRINFO(B,2)<N THEN ERROR(QSOLVE: B TOO SMALL);
	IF ARRINFO(X,2)<N THEN ERROR(QSOLVE: X TOO SMALL);
	BEGIN "QCOMP" REAL ARRAY LUA[1:N,1:N]; INTEGER ARRAY PS[1:N];
	DECOM(A,LUA,PS); SOLVE(LUA,B,X,PS);
	END "QCOMP"; END "QSOLVE";

INTERNAL PROCEDURE INVRT(REAL ARRAY A,AI);
		⊃ F-M matrix inversion (see p.78) without IMPROVE;
	BEGIN REAL N; N←ARRINFO(A,2);
	IF ARRINFO(A,4)≠N THEN ERROR(INVERT: [A] NOT SQUARE);
	IF ¬(ARRINFO(AI,2)=ARRINFO(AI,4)=N) THEN ERROR(INVERT: [AI] WRONG SIZE);
	BEGIN "ICOMP" REAL ARRAY LUA[1:N,1:N],B,X[1:N]; INTEGER ARRAY PS[1:N];
	INTEGER I,J; DECOM(A,LUA,PS);
	FOR J←STX N DO BEGIN FOR I←STX N DO B[I]←IF I=J THEN 1 ELSE 0;
		SOLVE(LUA,B,X,PS);
		FOR I←STX N DO AI[I,J]←X[I]; END;
	END "ICOMP"; END "INVRT";

INTERNAL REAL PROCEDURE DETERM(REAL ARRAY A); 
		⊃ Returns determinant of A using results of SOLVE;
	BEGIN INTEGER I,N; REAL DET; N←ARRINFO(A,2);
	IF ARRINFO(A,4)≠N THEN BEGIN OUTSTR(CRLF&"DETERM: [A] NOT SQUARE");RETURN(0) END;
	BEGIN "DCOMP" REAL ARRAY LUA[1:N,1:N]; INTEGER ARRAY PS[1:N];
	DECOM(A,LUA,PS); DET←1;
	FOR I←STX N DO DET←DET*LUA[I,I];
	END "DCOMP"; RETURN(DET) END "DETERM";

END "MATRIX";